Mon. Not. Ft. Astron. Soc. 000, [|-|l4| (2000) Printed 1 February 2008 (MN WF$l style file vl.4) 



Systematic uncertainties in gravitational lensing models: a 
semi-analytical study of PG1115+080 

HongSheng Zhao & Danny Pronk 

Sterrewacht Leiden, Niels Bohrweg 2, 2333 CA, Leiden, The Netherlands 



f^*^ ' Accepted Received ; in original form 

o ' 
o 

(N 

; ABSTRACT 

While the Hubble constant can be derived from observable time delays between images 

£/3 ' of lensed quasars, the result is often highly sensitive to assumptions and systematic 

uncertainties in the lensing model. Unlike most previous authors we put minimal 

restrictions on the radial profile of the lens and allow for non-elliptical lens potentials. 

£SJ . We explore these effects using a broad class of models with a lens potential i(j(rf(8)), 

which has an unrestricted radial profile but self-similar iso-potential contours defined 

by rf(6) — constant. For these potentials the lens equations can be solved semi- 

analytically. The axis ratio and position angle of the lens can be determined from the 

image positions of quadruple gravitational lensed systems directly, independent of the 

' radial profile. We give simple equations for estimating the power-law slope of the lens 

f^) , density directly from the image positions and for estimating the time delay ratios. Our 

method greatly simplifies the numerics for fitting observations and is fast in exploring 

the model parameter space. As an illustration we apply the model to PG1115+080. 

An entire one-parameter sequence of models fit the observations exactly. We show that 

q the measured image positions and time delays do not uniquely determine the Hubble 

5— i ' constant, 
-i— > ■ 

^ | Key words: gravitational lensing - quasars: individual: PG1115+080 - distance scale 

- dark matter - methods: analytical 
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C3 \ 1 INTRODUCTION 

Gravitationally lensed systems are powerful probes of galactic potentials and the scale of the universe. The advantage over 
the traditional stellar dynamical method is that from the image positions we can measure the shape and the mass of the 
dark halo well beyond the half-light- radius of a faraway lens galaxy, whether it is virialized or not. The time delay between 
two images is a measurement of the difference in the length of the two bent light paths, and scales with the distances to the 
lens and source. So the Hubble constant can be constrained once the redshift of the lens and the source and the time delay 
are measured. This way of getting Ho has the advantage that the underlying physics (general relativity) can be rigorously 
modeled. The limitation is that there is often a sequence of lens models that can fit the image positions. 

Presently about twenty strongly lensed systems are known, half of them being quadruple-imaged quasars and half being 
double-imaged quasars (e.g., Keeton & Kochanek 1996). We will concentrate on quadruple- imaged systems. They are better 
constrained than double-imaged systems, since the lens model needs to fit more image positions and also the ratios of time 
delays between any two pairs of images. Quadruple systems typically involve a quasar source well-aligned with the center of 
the lens potential well. Presently only two such systems (PG1115+080 and B1608+656) have accurately determined image 
positions and time delays. 

A significant amount of numerical computation is usually required to invert image positions to intrinsic parameters of the 
potential (cf. Schneider, Ehlers & Falco 1992). The degeneracy of the resulting potential is often not fully explored because 
of the need to cover a large parameter space, particularly for flattened potentials. Previous authors have often restricted their 
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studies to isothermal or power-law spherical models (Evans & Wilkinson 1998 and references therein), and elliptical models 
(Witt & Mao 1997 and references therein) and other simple models (Kassiola & Kovner 1993, 1995) with or without external 
shear. The fully general non-parametric method, e.g., the pixelated lens method of Williams & Saha (2000), is very powerful 
in demonstrating the complete range of the degeneracy in the lens models, but it involves significant amount of numerical 
computation and does not provide a clear insight to the relations between the characteristic parameters of the lens. For these 
reasons it is still desirable to find analytical, yet general potentials, which allow a quick exploration of the model parameter 
space. For example, it is of interest to generalize the analytical work of Witt & Mao (1997) to non-elliptical lenses. It would 
also be interesting to find analytical expressions for the time delay in these general lenses, which could help us to understand 
how the radial profile and lens shape affect the predictions on the Hubble constant. 

Here we study a broad class of analytical models with non-axisymmetric, non-elliptical shape and semi-power-law radial 
profile (§2), and show how to calculate the lens shape and radial profile parameters directly from the image positions (§3). 
We apply the models to PG1115+080 (§4) and show that the images can be fit perfectly by a large range of lens models (§5). 
We summarize our results in §6 and conclude with the implications on the Hubble constant. 



2 LENS EQUATION IN A GENERAL CLASS OF MODELS 
2.1 Decoupling of angular dependence and radial profile 

Any two-dimensional lens potential can be cast in the following form, 

ip(x, y) = i/>(u)), ui = u>(r, 9), (1) 

where tj) has the dimension of square arcsec, (x, y) defines a rectangular coordinate system (in units of arcsecs to the West 
and North of the lens galaxy center) and (r, 9) the corresponding polar coordinate system with 

(x, y) = (— r sin 9, r cos 9), (2) 

where 6 is the position angle, counterclockwise from North. Unless otherwise specified we shall follow the notations of Schneider 
et al. (1992). Here ui is defined to have the dimension of radius, so that constant u> curves correspond to equal-potential 
contours, and define the shape and flattening of the potential. The radial profile of the potential is a smooth function i/)(u>) 
of the radius u). We can also define 

m(uo) = w-j^- (3) 

as the mass (in units of square arcsec) enclosed inside the radius u> (in units of arcsec). For example, we have m = aip for a 

power- law model ip oc uo a with slope a. For a source at redshift z s and angular distance D os (z s ), and a lens at redshift zi and 
distance D i(zi), the physical mass M(u>) is related to m(w) by the scaling 

M(u>) = ttE c m(u>), (4) 
where 

s _ c Dpi Dos / 1 radian \ 2 ^ DoiD os ,^ 

c ~~ 4ttG D ls \ 206265 arcsec/ ~ D ls ^ ' 

is the critical density in units of Mq /arcsec 2 , and Di s (zi,z s ) is the relative distance of the lens and source and all distances 
are in units of parsec. 

A light ray from a source at (r s , 9 S ), being deflected to a direction (r, 9) by the lensing galaxy with potential ip(to) located 
at the origin, will experience a time delay At given by 

T 1 1 

At(r, 9) = h~ 1 T 100 (zi,z s ) -r 2 - rr s cos(6» - 9 S ) + -r 2 s - ip(v) , (6) 

where h is the Hubble constant Ho rescaled to 100 km/s/Mpc. A characteristic value for the time delay is 

noo(zi,z s )= Do ^ >osl + Zl , for Ho = 100 km/s/Mpc. (7) 

According to Fermat's principle, the images lie at the minimum of At, so the lens equation is given by 

-^-At = r - r s cos(6» - 9 S ) - m(u)d r lnw = (8) 
or 

-ItjA* = r s sin(0-0 a )-^^c> e lnw = O. (9) 
ro9 r 

Interestingly, the radial part m(u>) can be eliminated by simply combining the two equations, 
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r s sin(0 - S ) __ de inui , , 

r — r s cos(0 — 9 S ) rd r \nui' 

similar as in Witt & Mao (1997). This implies that there is a relation linking the image positions to the shape of the potential 
directly, independent of the radial profile. 



(16) 



2.2 Property of the image positions: the semi-hyperbolic curve 

For simplicity we shall concentrate on self-similar models with 

w(r,0)=r-/(0), (11) 

where f(9) defines the shape of the equal potential contours. In principle the shape function can be bi-symmetric or lopsided 
as long as the corresponding surface density is positive everywhere in the lens plane. To be specific, we will restrict our 
discussions to bi-symmetric potentials with an angular part 

/(0) = |l-<5cos20'|", (12) 

which is a three-parameter (/3,8,9 P ) function of the angle 0, where [3 is a constant, the parameter 8 is a flattening indicator, 
and 9 P is the position angle of a principal axis of the potential. The angle 9' is the azimuthal angle 9 except for a rotation 
with 

9' = 6-6 p , (x',y) = {r cos 9', r sin 8') (13) 

so that (x',y') defines a rectangular coordinate system with the axes coinciding with the principal axes of the lens. For 
example, the source at 

(x a , y s ) = (-r s sin0 s , r 3 cos0 s ) (14) 
in the original rectangular coordinate system would be at 

(x' s ,y's) = (r s cos 9' s ,r s sin6' s ), 8' 3 = 8 S - 9 P , (15) 

in the rotated rectangular coordinate system. 

With these we can compute the right-hand part of eq. [loj 

deinuj _ rfln/(0) 1 

rd r inu d9 v cot 9' + u tan 9' ' 

where 

A = 2P5, (17) 

are shape indicators just like (3 and 8. Substituting in eq. jj^, and rewriting the image radius rasa function of 9' , we find the 
images fall on a family of curves r = r(9') defined by 

r{0') = r s [cos(0' - 9' 3 ) + sin(0' - 9' s ) (y cot 9' + utan0')] . (18) 
An alternative expression for these curves can be obtained by expanding the sinusoidal terms in eq. [ll] so that 

r = (X a cos 9' + Y a sin 8') + ( -*L^ _ JL-) , (19) 

V cos 9' sin 9' J 

which is now linear in four new parameters (X a , Y a , X B , Y ); these parameters are related to the lens shape parameters (v, u) 
and the source position (x' s ,y' s ) by 

X a = (1 + v -u)x' s , Y a = (1 - u + v) y' s , X b = ux' s , Y b = vy' s . (20) 

These curves (cf. eq. |l8| ) have the nice property that they go through all image positions independent of the radial profile 
of the lensing galaxy. An example is the semi-hyperbolic curve in Fig. [I]. The curve is determined by the source parameters 
(r s , 0^), the lens shape parameters (v, u) and the lens position angle 9 P . The radial profile can take any general physical profile, 
isothermal or power-law. 

The boxiness parameter (3 is such that the shape function /(0) reduces to the usual elliptical form when f3 = | (Witt & 
Mao 1997). In the case that ip(uj) is linear in ui, the models reduce to the simple models of Kassiola & Kovner (1995) when 
[3 = 1. Interestingly, elliptical models with f3 = i have u — v = = jp = 1 (cf- ec i- |lj]), hence X a = Y a = 0, and eq. |l^ 
reduces to 

1=^-4 (21) 
x' y' 
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after applying x' — rcos8',y' — rsinff . This equation prescribes a hyperbolic curve, which is consistent with Witt (1996) 
and Witt & Mao's (1997) finding that all four image positions and the source position lie on a certain hyperbolic curve. A 
hyperbolic curve has a maximum of 5 free parameters, thus they cannot fit image positions of a general quadruple system; 
four image positions yield a minimal of eight constraints. Experience with fitting several quadruple systems (G2237+0305, 
CLASS1608+656, HST12531-2914) shows that models with \(5\ > 1/2 often give unphysical mass-radius relations m(w). We 
find that P = —\ gives a fair approximation to realistic models. These models have non-elliptical contours, and often yield 
physical density distributions. Fig. ^| shows that they also cover a sufficiently wide range of axis ratios for the potential and 
the density so that we can explore the shape of the lens galaxy in fitting the image positions. The expressions for the axis 
ratios of a power-law lens are given in Appendix A. 



3 RESULTS 

3.1 Lens shape directly from fitting image positions 

Our models can be used to fit image positions and derive lens shape parameters (6,/3) and the source position (r s , 9s), free 
from assumptions of the lens radial profile, but subject to the assumption that the lens' angular profile obeyes eq. n2j. The 
four unknowns can be derived from the four observed image positions, i.e., the eight observables (ri,9i) with i = 1,2,3,4. 
The position angle 9 P of the lens principal axis is treated as a free variable. 

The procedure is simple. First substitute the four observed image positions in eq. ^ to obtain the following four linear 
equations 

+ + -J— X b ?- 57 H = l, * = 1,2,3,4, (22) 

ri n n cos 9\ n sm 9\ 

of the four unknown parameters (X a ,Y a ,Xt,,Yb). After solving these, either analytically or numerically, these parameters are 
substituted in eqs. [IHJ and]!?] to yield the source position and the lens shape in terms of (r s , 9 S , 5,(3). In fact, we can recast 
eq. ^ to a set of four simple linear equations of four new unknowns (u, v,l/x' s , l/y' s ) by moving the terms x' s and y' s to the 
left hand side of the equations, i.e., 

X a /a;' s -v + u = l, Y a /y' s - v + u = 1, Xt,/x' s -u = 0, Yb/y' s -v = 0. (23) 
The lens shape parameters (S, (3) can then be computed from 

u + v 2(u — v) 

and the source position (r s ,9 a ) from 

r s = \J x' s 2 + y' a 2 , 9 S = arctan(2/s,a;' s ) + 8 P . (25) 

Thus we have effectively reduced the problem of fitting image positions to successively solving linear equations, which is a 
straightforward task. 

3.2 Non-parametric radial profile and power-law slope 

The radial part of the potential can also be extracted from eq. H without parameterization. At the positions of the images we 
have 

m(oj) = r 2 — rr s cos(# — 9 S ) — x 2 + y 2 — xx s — yy s , u> = rf(9), (26) 

where we have used d r Into — r _1 . Thus we have obtained a mass-radius relation directly from the observed image positions, 
assuming that the source position (x 3 ,y s ) and the flattening and position angle of the potential (S,9 P ) have been determined 
by fitting the curve (cf. eq. |l8| ). Interestingly, in the limit that the source is at the center of a circular lens, we have r s — > 0, 
,9^0 and m — > r 2 . 

It is useful to characterize the radial profile of a lens galaxy, which is generally not a power-law, by an effective power-law 
slope which varies with the radius to except for scale-free models. There are several ways of estimating the characteristic 

power-law slope. Taking any two images i and j, we can form a characteristic power-law slope ay from the mass m; and rrij 
at the image radii u>i and u)j with 



21og^ + log 



log mi— logm, ° T i ' ° l-r 3 r/ r cos(8 i -8 a ) 



(27) 



13 " log ^+ Hog \zlZ£~-e') ' 
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where we have used the mass-radius eq. [26| 

Alternatively we can estimate the power-law slope from the observed time delay. First we rewrite the time delay eq. ^| 
so that the time delay, ty , between any two images i and j is given by 

— = ijfci + yf) -\(x* + y)) - (xi - Xj)x a - (yi - yj)y s - (ipi - tpj), (28) 

TlOO 2 1 

where ipi — iftj is the difference in the lens potential between the two images. Rewriting eq. ^ we form a new estimator a*j 
for the power-law slope with 



0? + Vi) - (xj + y)) - (Xi - Xj)x s - ( yi - yj)y a 



(29) 



T ioo 



V iH-ti l( a; 2 +y | ) _l (x 2 +y 2 ) _ (x ._ x . )a . a _ (y ._ y . )ya _M 

where we have used eq. ^(j. 

Thus we have two direct estimators ay and aL of the radial profile, computed from the observed images and time delays. 
In the limit of scale-free power-law models ay = ay = oi(lu) = cat. So the deviation from scale-freeness can be estimated by 
taking the differences such as ai2 — «34, a*2 — 034, or ai4 — 0*4. 



3.3 Hubble constant and time delay ratios 

To apply the above estimates of the power-law slopes, we have assumed that we know the rescaled Hubble constant h from 
independent observations. Alternatively the Hubble constant Hq = lOO/i can be estimated from the time delay ty between 
two images. Letting ah = ay, we get 



The time delay ratio can also be predicted with 

t .. (3 - ^7) [(*? + Vi) - + Vj)] ~ (l ~ ^7) - + (Vi ~ Vo)Vs] ■ 



which is obtained by rewriting eq. H. 



(31) 



4 APPLICATION: THE SURFACE DENSITY AND TIME DELAY MODELS FOR PG1115+080 
4.1 Data 

As a simple application, we model the image positions and time delays of the well-studied quadruple system PG1115+080. 
This system has been extensively studied ever since the first models by Young and collaborators (1981), and has received 
closer attention after Schechter et al.'s (1997) measurements of its time delay. All models, except those of Saha & Williams 
(1997), adopt elliptical/circular shapes and a few common radial profiles, with the models of Keeton & Kochanek (1997) 
and Impey et al. (1998) being the most comprehensive. Only one year after the discovery of the legendary double-image 
radio-loud quasar Q0957+561, this system was identified as a multiple-imaged system by Weymann et al. (1980) in their 
survey of nearby bright QSOs. It is now known to consist of four images with the names A\, A2, B and C (with flux ratios 
about 4 : 2.5 : 0.7 : 1, cf. the HST observations of Kristian et al. 1993) of a radio-quiet QSO at redshift z s = 1.722. The 
images Ai and A2 are within 0.5"; see the inset of Fig. |l| Interestingly, the lens galaxy is also one of the bright members 
of a galaxy group (N ~ 10) at redshift zi = 0.310, first mapped by Young et al. (1981). The center of the group is to the 
south-west of the lens, roughly at r g = (20" ± 2") and g = (-117° ± 3°). The lens galaxy has been resolved by both HST 
and the 8.2-m Subaru telescope in 0.3" seeing. It appears to be an early type galaxy with a de Vaucouleurs profile and a 
half-light radius of 0.55". There is no sign of differential dust-extinction in the lens galaxy. While NICMOS observations by 
Impey et al. (1998) show no flattening for the lens, ground infrared images by Iwamuro et al. (2000) found it to be an El 
galaxy elongated towards 9 P ~ 65°. Both observations reveal a 1" infrared Einstein ring connecting the four images, which 
is thought to be the infrared image of the QSO host galaxy. PG1115+080 is also one of the two quadruple systems where 
the time delay between images has been measured, the other one being the radio-loud quasar B1608+656 from the CLASS 
survey (cf. Fassnacht et al. 1999). Although two different sets of values are quoted in the literature (Schechter et al. 1997, 
Barkana 1997), the leading image is the furtherest image (the image C), and the innermost image (image B) arrives last. The 
time delay ratio tabc = tAc/tBA, for the delay between image C and Ai + A2 vs. image B and Ai + A2, provides an extra 
discriminator of the models; the images Ai and A2 are within 0.5" of each other, and the small relative delay is undetected. 
Schechter et al. first reported tabc = 0.7 ± 0.3 from their photometric monitoring program in 1995-1996. Later analysis by 
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Barkana (1997) found tabc = 1.13 ±0.18, after taking into account correlations of errors amongst the time delays. The delay 
between images B and C, tsc ~ 25.0 ± 1.7 days. 

Here we illustrate the application of our models to the most recent data from Impey et al. (1998) of PG1115+080. We 
do not attempt to model B1608+656 because the lens appears to be a merging pair of galaxies, and the morphology is too 
complex for our model. We denote with the index i = 1, 2, 3, 4 the four images Ai, A2, B and C. All results are quoted for a 
standard flat universe without a cosmological constant (f2, A) = (1,0). Table [l] gives the relevant quantities to calibrate our 
results to other universes. The angular distance from redshift z\ to zi in a universe of a matter and vacuum density (f2, A) 
times the closure density is generally given by 

d(zi,2 2 ) = — ; — c _ sinh \VsZ [ dz\A + n k (i + z) 2 + (i + z) :i nr^\, n fc = i-A-n. (32) 

Ho(l + z 2 )V^c L h x i 

The predicted Hubble constant should be reduced from the value achieved with the standard (Q, A) = (1,0) universe by 3% 
in the presently favored A-dominated universe with (f2, A) = (0.3, 0.7) from surveys of distant supernovae. 



4.2 The source position, the lens shape and mass inside images 

First we solve for the lens shape and source position from the linear equations ^ and ^3| The solutions for PG1115+080 are 
given in Table |^, sorted according to the value of (3 or 8 P . The resulting potential model with (3 — — | or 8 P = 60.5° has a 
flattening of between E0 and El, and interestingly the lens principal axis points towards the location of the galaxy group in 
the lens plane. Models with other values of /3 or 9 P will be discussed in section 5. 

To proceed with determining the lens mass at each image position, we substitute the now known flattening parameters 
(5, 9 P ) and the source positions (r Sy 8 s ) in eq. |26| to predict four independent data points (oji,mi) with i = 1,2,3,4 in the 
radius vs. the enclosed mass plane. Fig. |^ shows the predicted lens mass enclosed at the four image positions. Note that the 
mass rises faster than the light, implying a growing dark mass component at large radius; the light distribution is modeled as 
an observed de Vaucouleurs r*-law with a half-light radius of 0.55". 



4.3 Piecewise power-law model 

So far we did not enforce any strict parameterization of the radial profile. We only restrict the profile to be of the form of 
eq. |ll]; in practice, this is an insignificant restriction. In the following sections we show several ways of modeling the radial 
profile assuming an isolated lens. None of the models is completely satisfactory. 

First we use a minimal model, which assumes that the mass-radius relation is a piecewise power-law, that is, we connect 
a straight line through two images in the log m vs. log uj plane. The piecewise values for the power-law slope and axis ratios 
are given in Table ^. 

The Hubble constant Hq — 100ft can be estimated by normalizing the time delay to the observed value tAC (Schechter 
et al. 1997), 



, T 100 

h ■ 



'it AC 

where 



-/',,-- (1- — J.S'n 
ai4 



(33) 



Pij= [(*?+»?)-(*? + «?)] (34) 
depends on the image positions (xi,yi) alone and 

Sij = 2Pij - 2 \(xi - Xj)x s + {m - yj)y s ] (35) 

depends on the source position (x B ,y s ) as well. The model yields a Hubble constant Hq ~ 30 km/s/Mpc, much lower than 
determined by other authors (e.g., Impey et al. 1998). The low Ho is a result of the high value for the power-law slope 
cxac = 1-6. Other values are given in Table ^ for various image pairs and observed time delay. The time delay ratio can also 
be predicted with (cf. eq. ^]J) 

_ tu -J'u+Cl-^Su 

TABC = " = 7 \ — \ ■ (36) 

*32 -P 32 + (1 - -L) 5 32 

Substituting the slopes oac and oiba to eq. ^ we find the ratio between images tabc = tAc/tBA ~ 0.65. 

Note that the time delay predictions here are robust and independent of details of the density, since the discontinuity in 
the density is completely smoothed out in the lensing potential. 
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4.4 Single power-law model 

The piecewise-power-law model above necessarily implies a discontinuous density profile. This could be cured if we enforce 
a single power-law model, that is, we fit a straight line to the four points in the log m- log u plane. This would give us a 
power-law slope a — 1.38. Fig. |l| shows our model surface density contours. Similar to the non-parametric models of Saha & 
Williams (1997) we find a peanut-shaped lens. 

We can estimate the goodness of the fit by recomputing the image positions from a given mass model. The general 
procedure of simulating images of our theoretical lens model is as follows: First combine eq. ^ with the power-law radial 
profile to get 

r 2 - rr s cos(6 - 9 S ) = m{uj) = b r > w = r 1 1 - 5 cos(2<9 - 29 p ) f . (37) 

Then upon substitution of eq. ^ to eliminate r, we obtain a one-dimensional non-linear equation for the image position angle 
9, which can be solved easily numerically. As a comparison, one would be dealing with a minimum-finding or a root-finding 
numerical problem in a two-dimensional plane (r, 9) in the general case without the separation of the angular vs. radial part. 
For our model of PG1115+080, the image solutions are shown in Fig. ^ as dashed circles, together with the input observed 
image positions (solid circles). The predicted four images are off by about 60 milli arcsecs, a residual which is inconsistent 
with the ~ 3 milli arcsecs accuracy of Impey et al. positions, and is marginally consistent with earlier data by Kristian et 
al. with a quoted error of 50 milli arcsec for the lens galaxy and 5 milli arcsec for the images. The images Ai and A2 are at 
two sides of the critical curve (the line of infinite amplification in the source plane), hence are highly amplified with opposite 
parity. The time delay ratio can be estimated with eq. BU assuming a constant power- law slope a — cxbc = 1.38. This yields 
tabc ~ 1-3, close to the Barkana (1997) value of 1.13 ± 0.18. The large difference here is due to the large residual in terms 
of fitting the images Ai and A2 with a straight power-law. 

We also compute the amplification patterns by taking double derivatives of the time delay surface. The circles in Fig. [j] 
show the observed (solid circles) and predicted (dashed circles) fluxes of each image and the source (half-closed circle), with 
the area of each circle in proportion to the flux. The B to C ratio is well-reproduced and the predictions for images Ai and 
A2 are also consistent with observations at about the 0.3 magnitude level. 

Finally, if we put a host galaxy around the QSO, the model predicts that the image of the host galaxy will be stretched 
into an arc. We see a nearly closed ring. This agrees very well with the diffuse ring that Impey et al. discovered in their 
NICMOS images. The ring maps back to the source plane as a disk with an area of ~ 0.03 square arcsec. 



4.5 Double power-law model 

Alternatively we can fit a smooth, five-parameter lens model with a lens potential 

ui = r |1 - <5cos(26» - 28 p )\ 13 . (38) 

The corresponding mass profile is given in Appendix B. This model assumes the lens potential (as well as the lens mass) 
increases like a double-power- law with an inner slope a in and outer slope a out . The transition is defined by the normalization 
constant Co, the radius r — ao and the sharpness parameter n; bigger n corresponds to sharper transition. When fitting the 
mass model to the four points in Fig. ^, it turns out that the inner slope a4 n is fully unconstrained. The other four free 
parameters (co, ao, ot out , n) are determined from fitting the four data points; the procedure is explained in Appendix B. The 
values of (co, ao, a ou t) ~ (0.6, 1.1, 1.63), nearly independent of the value for the inner slope am- Note a out ~ «ac, consistent 
with the fact that the mass at Ai, A2 and C follows nearly a power-law (cf. Fig. ^|). The value of n increases from ~ 15 to 
~ 39 for ai n = — 2. That n ^> 2 means that the density changes sharply at the transition. All fits have zero residual and 
predict nearly identical mass profiles at radii between images B and C. They differ only in the mass profile inside image B, 
where we have no direct constraint on the dark matter profile. For all purposes it is sufficient to set n = 20 so that aj n ~ 2, 
in which case the model has a finite core at small radii. From the potential model we can compute the dimensionless surface 
density contours (cf. Fig. Q). Near the position of the images, the density has a flattening of E2-E3, flatter than the potential, 
as expected. The potential model can also be substituted in eq. |^ to predict the time delay contours, shown in Fig. ^. The 
observed images C, A2, Ai, B (the solid circles) are exactly the valley, peak, saddle, valley points of the delay contour, where 
the theoretical images should lie according to Fermat's principle. The images Ai and A2 have nearly the same arrival time. 
Substituting the lens potential in eq. ^S] we can predict the time delay ratio for the double-power-law model. The result, 
tabc ~ 0.65, is nearly the same as the piece-wise power-law model. These predicted ratios are in good agreement with the 
earlier measurement of 0.7 ± 0.3 of Schechter et al. (1997). 

We can recompute the image positions from the double-power-law model by solving 



uj = r |1 - Scos(29 - 20 p )f , (39) 
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and eq. [L8|. 

Unlike the single-power-law model it predicts five images (cf. Fig. ^J). Four predicted images fall exactly on the observed 
positions. But the model predicts an extra image near image B. This can be identified with the usually highly demagnified 
fifth image that is sometimes observed in lensed systems. The extra image the model predicts is, however, too bright to be 
consistent with observations. This is a direct consequence of the fact that k at its radius is nearly constant and ~ 1. The 
amplification can be calculated from k as (1 - k)~ 2 ~ 25. 

All images lie on the curve defined by eq. The extra image arises because the model density profile (cf. Fig. ^ is 
non-monotonic near 1" radius, the transition radius ao of the double-power-law potential model. A wiggle in density happens 
when the transition parameter n > 2, i.e., a sharp transition of the potential. 



5 DISCUSSION 

5.1 Power-law slope of PG1115+080 vs. Ho and time delay ratio 

The previous sections have illustrated the procedure for fitting known quadruple image systems with our models, and have 
also revealed some puzzling problems with PG1115+080, for example, the low and non-unique Ho and the peculiar and 
non-unique lens density. Fig. ^| shows the values for a predicted from either the mass-radius relation (cf. eq. |2t| ) or from the 
time delay (cf. eq. ^). Two comments are in order: (1) There is a large spread with the power-law slope a depending on how 
it is predicted. To fit the observed time delay within 2a with a reasonable Ho (between 50 and 100 km/s/Mpc), we obtain 
only a loose constraint 0.3 < a < 1.6. (2) The rise of the power-law slope a with radius in the piecewise-power-law model 
is somewhat unphysical; realistic models typically have a monotonically steeper density profile (smaller a) with increasing 
radius. These results are consistent with the finding of a wiggle near the transition radius (1") in the surface density of the 
smooth double-power-law fits. Nevertheless the density is positive everywhere. 

A detailed treatment of these problems should include the effect of the neighboring group. Such models are presented 
elsewhere, since they are beyond the interest of this paper on methods for quadruple systems in general. Nevertheless some 
further study of the model parameter space is clearly needed. Here it is sufficient to comment on the range of the predicted 
model parameters, in particular, the power-law slope and the value of Ho when we vary the parameter f3 or 8 P . 

Fig. |7| sh ows the parameter space of the isolated lens model by varying the orientation of the lens principal axis 9 P (cf. 
eq. ^ and |l3|) . It turns out a one-parameter family of models fits the observations; we effectively vary j3. In fact, 9 P ~ 61° +4/3, 
and A = 2/3S ~ § to a fair accuracy. We find that models with /3 < —1/3 will generally result in negative density models, 
models with j3 > have a radially-increasing density, and models with |/3| — ► also produce a nearly constant axisymmetric 
density, hence implying an unphysically small Hubble constant. Only in the range — j < (3 < do we find physical models 
with a reasonable flattening. The source positions for different models are shown in Fig. 0. Table ^ compares the predicted 
parameters of the j3 = — ^ or 9 P — 60° model with the f3 = — | or 9 P = 60.5° model. Fig. pTshows the predicted radial profile 
for both models. Both models would have a significant residual if we fit the image positions with a straight power-law, or 
strange wiggles and extra images if we fit a double-power-law. We find this is generally the case with all isolated lens models 
without the shear from the neighbouring group. 

Fig. |^ shows that there is significant degeneracy with the value for Ho, depending on the adopted model and the adopted 
time delay. The most reliable estimate is from the delay between the innermost image B and the outermost image C using 



, T 100 

h ■ 



2ti 



-Pli + ( 1 - —) S U ~ ^32 + f 1 - —) S 32 

\ Ql4/ V Q32 / 



(40) 



where we neglected the time delay between the images A\ and Ai- The Ho thus predicted scales approximately as 60(2 — 
a) km/s/Mpc (cf. Fig. []), where a = oibc is the average power-law slope. The factor 2 — a is the average power-law slope of 
the model surface density, and is not well-constrained at the radius of the images. Models with f3 ~ — | or 6 P = 60° predict 
a steeper density profile, and closer to isothermal than models with j3 ~ — | or 6 P — 60.5°, hence yield a more plausible Ho 
around 60 km/s/Mpc. Models with very shallow profiles (a > 1) are unfavored by the consensus value for Ho- 

In all our isolated lens models, the delay ratio is closer to the Schechter value, while the Barkana value seems to be the 
widely accepted one. This explains why we see in Fig. |^ a tighter prediction of Ho using tAC and tsA of Schechter values 
than the Barkana values. Interestingly in the limit that the isothermal model is applicable, a = 1, the time delay ratio can 
be estimated from the image positions directly (i.e. without computing the source position) with the following equation: 

P r 2 — r 2 

f abc = tt— = ~k % ~ 0.65, for isothermal lens, (41) 

P32 r 2 - r 2 

where rj is the observed radius of the image i from the lens center. The ratio increases when we introduce the shear from the 
nearby group. 
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5.2 Rotation curve and mass-to-light ratio 

Fig. |l(j compares the predicted velocity dispersion for the lens with the observed value. The dispersion is predicted with the 
following formula, 

G 2 -V 2 — V 2 — Dos ^ X lo5km / s ) 2 1 radian 
1 2rt ' Di 3 2tv radian 206265 arcsec ' 

where the predictions are made from each image i = 1,2,3,4 (strictly speaking the formula is valid only for a singular 
isothermal lens). Applying this for each value of the lens position angle 9 P we get the full range of possible values for the 
dispersion. The difference of the dispersion among the four images shows the deviation from an isothermal model. On average 
the predicted lens galaxy velocity dispersion, 200-300 km/s, suggests that the lens is close to a massive L* galaxy. The 
observed dispersion a = (281 ± 50) km/s (95% confidence, shaded region) from the Keck spectrum of Tonry (1998). This 
value is comparable to the theoretical models, but on the high side. A similar problem has been noted in previous models 
(Schechter et al. 1997). Tonry noted that the discrepancy might be due to a steep radial fall-off of the velocity dispersion or a 
simple over-estimation of the observed dispersion; the spectrum of the faint lens was made from subtracting two Keck spectra 
taken with a l"-slit, one passing the lens and the image B, one passing the images A\ and Ai. 

Fig. |ll] compares the increase of the lens mass from the innermost image B to the outermost image C vs the increase of 
the lens light. The mass ratio is predicted using 

M£. = T± = ( r _±\ asc / 43 ) 
M B m 3 \r 3 7 

where the power-law slope cxbc depends the principal axis 9 P . The light ratio 



Lb 



J Q 4 dr2-Kr exp 


-7.67 


J Q 3 dr2-Kr exp 


-7.67 (i) 4 



= 0.55", (44) 



where the observed de Vaucouleurs law is used for the projected light. The fact that -j^- > > 1 implies that the mass 
grows faster than the light at these radii for all these models, consistent with a large amount of dark matter between 0.8 and 
1.5 arcsec. 



6 CONCLUSION 

Here we itemize our main results. 

(i) We found a very general class of lens models that allow for non-elliptical and non-scale-free lenses. These models 
include previous isothermal elliptical models as special case. We can derive the radial mass distribution of the lens in a 
non-parametric way. We can study the deviation from the usually assumed straight power-law profile. We also give simple 
formulas for computing the time delay ratios and for estimating Hq. 

(ii) The models are very easy to compute and can be used for efficient exploration of a large parameter space because 
the lens equations can be reduced to a set of linear equations. 

(iii) We have applied the models to PG1115+080, and have explored a large parameter space of isolated lens models. All 
models using piece- wise-power-law or double-power-law fit the positions of the images exactly. The models can also reproduce 
the flux ratios between the images and the stellar velocity dispersion of the lens approximately. 

(iv) Our models are consistent with a dark halo up to a radius of three times the half-light-radius of the lens. The enclosed 
mass increases much faster than the enclosed light as we move radially from the innermost image to the outermost image. 

(v) We reconfirm earlier results by (e.g. Schechter et al. 1997) that the principal axis of the lens potential points, within 
a few degrees, to the external group, and is consistent with the observed value (Iwamuro et al. 2000). 

(vi) Our models do not yield a unique prediction of Ho, because it is sensitive to the lens power-law slope a (cf. Fig. |^), 
which appears to have a large spread, depending on how it is predicted. The power-law slope is also sensitive to small 
uncertainties of the position angle of the lens 9 P . A change of 9 P by a few degrees from the observed value can change a from 
to 2 and H from 100 km/s/Mpc to 0. 

(vii) Our models predict consistently a low time delay ratio tabc, which fits only the Schechter value. This and the 
peculiar oscillation in the predicted density profiles, we believe, are due to the neglect of the external group. This will be dealt 
with in detail in a follow-up paper. 

In summary, our semi-analytical lens models allow us to explore a large range of physical lens density distributions. We 
find that there is a large systematic uncertainty in the lens models from fitting image positions and time delays of PG1115+080 
and isolated lens models are always unsatisfactory for this quadruple system. 

The authors thank Tim de Zeeuw for a careful reading of the manuscript and the referee for a constructive report. HSZ 
thanks Paul Schechter for discussions and encouragement. 
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Table 1. Effects of different cosmology on physical distance, velocity and time delay scales 



(n,A) 

kpc/l" h 

V 2 (km s- 1 ) 2 /!' 



(1.0 , ) (0.4 , 0.6 ) 

2.808 3.127 

316 2 305 2 

31.98 33.26 



(0.3 , 0.7) 



3.195 
303 2 
33.37 



(0.2 , 0.8) 



3.271 
299 2 
33.38 



T 100 (days per sq. arcsec) 



Table 2. Results for the models with the boxiness parameter f3 = 



1 and = -| 



Parameters /3 = — g 

PA 9 P 60.5° 

{r 3 ,0 s ) (0.09", 35°) 

qj, 0.92 

q K 0.73-0.77 



60.0° 
(0.03", 33°) 
0.86 
0.20-0.60 
(0.51,0.99,1.31) 




(a fl A,«BC.«ylc) (1.1,1.4,1.6) 

M/L/h 27.0 

tabc 0.6 

H 20-50 



26.3 
0.6 
50-90 
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Figure 1. shows the relative locations of the four QSO images and the lens galaxy (filled square at the center) in the PG1115+080 
system. The inset shows the relative locations of PG1115+080 and the neighbouring galaxies (Gl to Gil in the notation of Impey et al. 
1998). The equation [is] for self-similar models also predicts a semi-hyperbolic curve passing through the observed images (solid circles), 
the predicted images (dashed circles) and the source (smaller half-closed circle) and the lens (the filled square) with asymptotic lines 
always being parallel to the symmetry axes of the potential. The area of the circles is in proportion to the flux. Overplotted are the 
predicted surface density contour maps (solid contours for density being 1 and 2 times the critical density) and critical curves (dashed 
contours) in the single-power-law model with (3 = — i. The nearly closed ring is predicted for a uniform disk around the QSO. 
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Figure 3. The predicted lens mass (in units of lO 11 ^ 1 Mq) inside the radii of the four images. The predictions depend on the shape 
parameters of the lens, and are made here for = — g (principal axis e p = 60.5°) and /3 = - \ (6 P = 60°). Also shown is the observed light 

profile (arbitrary flux normalization) with a de Vaucouleurs r~ 4-law and a half-light radius of 0.55". Note a radially-rising mass-to-light 
ratio. 
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Figure 4. Similar to Fig. |l| but for the double-power-law model with f} = — |. Overplotted is a solid contour, showing the predicted 
surface density at the critical density. Also shown is the predicted time delay "finger print" for the double-power-law model with j3 = — g . 
The contour passing close to the separatrix (image B) is split to two disconnected contours. The extra fifth image (the isolated dashed 
circle near the image B) is a result of the non-monotonic model density, and it sits on a local extreme of the delay contour. The other 
contour is a contour of constant time delay ratio txc/^BX = 0.7 of any point X and images B, C. 
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APPENDIX A: POWER-LAW MODELS 



Within a certain range of radii, realistic lens mass distributions can be approximated by the following bi-symmetric power-law 
of the radius to 



m(w) 



dip 



w = rf(6) = r\l - 5cos(28- 26 p )f , r 



< a < 2, 



(Al) 



where bo is the characteristic deflection strength at the characteristic length scale ro (one arcsec) of the lens system, and 6 P 
defines the principle axis of the lens. The corresponding lens potential is related to the enclosed mass m(u>) by 



tp(oj) = const + 



m(w) 



The axis ratio (or its inverse) of the bi-symmetric potential 

q ^=(L=iy *>1-A, A = 2/35, 
and the axis ratio of the density 

V a 



\1+SJ 



2-Q 



1)8 



1- A(1 + 2q" 1 ), 



where the approximations are valid if the flattening is small. The surface density k of the model is given by 
K = ^br a ~ 2 g a0 (ag 2 + ([3 2 a - (3)g 2 + /3g"g), 
where g = 1 - <5cos(26> - 29 p ), g' = 2(5sin(26» - 26 p 



-45 2 cos{29- 29 p ) 



(A2) 



(A3) 



(A4) 



(A5) 



APPENDIX B: DOUBLE-POWER-LAW MODELS 

For the double power-law lens potential 



if)(ui) = CO 



\ao/ L \ ao 



the corresponding mass profile is given by 
dtp(u;) 

m(uj) = LU — p — - — a{UJ)1p(LO), 

duo 



where 



a(u>) 



(Bl) 



(B2) 



(B3) 



V ao / 

Generally speaking, a double-power-law fit (cf. eq. ^) involves searching for solutions of the following equations in a 



five-parameter space (n, £o, ce ou t, 7o, ao): 
log mi = [f + (1 + 7o)« otl tloga; l ] - ^1 + 
where 



70«out 



log 



1 + 



(-)' 

\aoJ 



log 



■ 7o ' 



(- 

Vao 



i= 1,2,3,4, 



7o = 1, Co = log c + log a,out — (1 + 7o)a otl i log ao, 



(B4) 



(B5) 



and {uJi, mi) for i — 1, 2, 3, 4 are the four data points in the mass-radius diagram (cf. Fig. ^). Note that the equations B4 are 
linear in (£,o,a out ), so these two variables can be eliminated by Gaussian substitution. If we fix n, say n = 20, then we are 
left with two equations for two variables (70,00), which can be solved with a two-dimensional iterative root finding routine. 
The solution converges rapidly from an initial guess of 70 = and ao = 1". 



This paper has been produced using the Royal Astronomical Society/Blackwell Science MpX style file. 
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Figure 5. The predicted run of the dimensionless density k along (solid) and perpendicular to (dashed) the principal axis 8 P = 60.5° 
(i.e. /3 = — A) for the double-power-law model. The model is positive everywhere, but has a peculiar wiggle at about 1". 
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0.5 1 1.5 2 



Figure 6. The predicted power-law slope of the lens galaxy, otij from fitting the mass-radius relation (the horizontal axis, cf. eq. |27]), 
and a*^ from fitting the time delay (the vertical axis, cf. eq. ps| ), where the ±2<r range of Barkana time delay and 50 < Hq < fOO are 
used with higher Ho corresponding to lower ct\-. The predictions are made for fitting the pairs (Ai, C), (B, C), and (A2, B). 
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Figure 7. The parameter space of the isolated lens model. The four panels show (in clockwise direction) the model boxiness parameter 
f3, the axis ratio <jy, of the potential and q K of the density (cf. eq. |As| and eq. |A4| ), the power-law slope a (cf. eq. p?] ), and the delay ratio 
y ABC (cf- e 1- The shaded bar at the bottom of each panel indicates the range of physically allowed models. 
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Figure 8. the predicted source position for the sequence of models with the principal axis in the range 58° — 62° (from the upper left 
to the lower right), which corresponds to j3 = —0.8 to /3 = 0.25. The source lies on a nearly straight line passing the origin, where the 
lens is. 
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Figure 9. The dependence of Hq on the power-law slope a, where a is determined by the ratio of increase in mass and increase in 
distance. The short vertical line to the left shows the slope of increase in light log(Lc /Lb)/ log(r£ /rg). The curves show the predictions 
for Ho using *ac> tgci and tgA measured by Schechter et al. (1997), solid line, and Barkana (1997), dashed line. The shaded region 
shows the consensus value for Hq. 
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Figure 10. The predicted velocity dispersion of dark matter in the lens galaxy for a series of lens models specified by different lens 
position angle 8 P . The values are calculated from the image positions of A\ (solid circle), A2 (open circle), B (triangle) and C (starry 
symbol) respectively. The shaded region shows the 95% confidence range from the observed average dispersion of stars in the lens (Tonry 
1998). 
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Figure 11. Similar to but for the predicted ratio of the lens mass inside the outermost image C and the innermost image B, 
Mc/Mg. Note that the ratio is always higher than the ratio of the light inside the two images, Lq/Lb (the dashed line); the latter is 
calculated from the observed de Vaucouleurs law. 
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